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Abstract 

Monte Carlo simulations within the grand canonical ensemble are used to obtain 
the joint distribution of density and energy fluctuations pl(p,u) for two model fluids: a 
decorated lattice gas and a polymer system. In the near critical region the form of pl(p, u ) 
is analysed using a mixed field finite-size-scaling theory that takes account of liquid- 
vapour asymmetry. Field mixing transformations are performed that map pl(p,u) onto 
the joint distribution of critical scaling operators p*(x,y) appropriate to the Ising fixed 
point. Carrying out this procedure permits a very accurate determination of the critical 
point parameters. By forming various projections oip k {x, y), the full universal finite-size 
spectrum of the critical density and energy distributions of fluids is also obtained. In the 
sub-critical coexistence region, an examination is made of the influence of field mixing 
on the asymmetry of the density distribution. 

PACS numbers 64.70F, 05.70.Jk 



1 Introduction 



Over the years, finite-size-scaling (FSS) techniques have proved themselves an indispensable 
tool for computer simulation investigations of critical phenomena in model systems, facilitating 
accurate estimates of infinite volume quantities from simulations of finite size. Of the many 
previous FSS simulation studies that have been performed [0J, most have focussed on critical 
phenomena in lattice-based magnetic spin systems such as the Ising [0, |3|, 4 ||, XY || and 
Heisenberg models || [7) . Among the specific approaches employed to study such systems, use 
of the order parameter distribution function has proved itself one of the most powerful. The 
FSS properties of the order parameter distribution are now routinely employed in studies of 
magnetic systems, facilitating both the accurate location of the critical point and the detailed 
elucidation of its character || . 

Only comparatively recently has attention turned to the task of applying FSS techniques 
to the simulation of critical fluids. Work to date has concentrated on attempting to carry 
over to fluids the order parameter block distribution techniques developed in the magnetic 



context [g, [10], |TTJ. In the process, however, it has been found necessary to generalise the 
FSS equations to take account of the absence in fluids of the energetic ('particle-hole' fT2|l ) 
symmetry that prevails in magnetic systems such as the Ising model. Although this reduced 
symmetry is believed to have no bearing on the universal critical behaviour of fluids, (which for 
system with short-ranged interactions correspond to the Ising universality class), it has long 
been recognised that it leads to certain non-universal effects. The most celebrated of these is 
a weak energy-like singularity in the coexistence diameter on the approach to criticality, the 



existence of which constitutes a failure for the law of rectilinear diameter |T3| 



Within the renormalisation group framework, the critical behaviour of a given system is 
characterised by the values of its relevant scaling fields specifying the location of the effective 
hamiltonian with respect to the fixed point |IJj]. In models of the Ising symmetry, these scaling 
fields are simply identifiable with the thermodynamic fields, namely the (reduced) temperature 
and the applied field. By contrast for fluids, the absence of particle-hole symmetry implies 
that the scaling fields comprise mixtures (i.e. linear combinations) of the temperature and 
chemical potential. As a consequence of this 'field mixing', the fluid scaling operators (the 
quantities conjugate to the two relevant scaling fields) are also predicted to differ from those 
of the symmetric Ising systems. In the Ising model, the scaling operators are simply the order 
parameter (i.e. magnetisation) and the energy density. In fluids however, they are expected 
to be linear combinations of the order parameter (particle density) and the energy density. 

The modified forms of the scaling variables have recently been incorporated within a FSS 
theory describing the interplay of energy and density fluctuations in near-critical fluids [|lfj, [LI]]. 
This theory provides a potentially powerful framework for the detailed simulation study of 
critical phenomena in fluids. Hitherto, however, only a limited appraisal of the theory has 
been performed. Attention has focused on one predicted manifestation of field mixing, namely 
a finite size correction to the limiting form of the critical order parameter distribution. The 
existence of this correction has indeed been confirmed in detailed Monte-Carlo simulation 
studies of both the 2D Lennard- Jones fluid fll| and a 2D decorated lattice gas model jlGfl . To 
date, however, the full extent of the claims embodied in the mixed field FSS theory have not 
been closely scrutinised. 

In the present paper we address this matter with simulation studies of two critical fluid 
systems: a decorated lattice gas model and a polymer model. The paper is broadly organised 
as follows. We begin by providing a short resume of the mixed field FSS theory for the near- 
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critical density and energy fluctuations of fluids. We then present simulation measurements 
of the joint distribution of the density and energy fluctuations pl{p,u) at the liquid-vapour 
critical point of a 3D decorated lattice gas model. Field mixing transformations are performed 
that map pi(p,u) onto the fixed point distribution of scaling operators p*(x,y) appropriate 
to the Ising universality class. Effecting this data collapse yields estimates of the field mixing 
parameters that control the degree of field mixing, the values of which are found to be in 
excellent agreement with analytic calculations. Application of field mixing transformation to 
p*(x,y) are also used to generate the full universal finite-size spectrum of the density and 
energy density distributions of fluids. This analysis reveals, in particular, that (compared to 
models of the Ising symmetry) the presence of field mixing radically alters the limiting (large 
L) form of the critical energy distribution. 

Consideration is then given to the role of field mixing in the sub-critical two phase regime of 
the decorated lattice gas model. An examination is made of the effect of applying field mixing 
transformations to the coexistence density and energy density distributions. For sub-critical 
temperatures down to approximately 0.9T C , it is found that the observed asymmetries of the 
coexistence density distributions are well accounted for by the linear mixing of the energy 
density into the ordering operator. 

Finally we apply the techniques developed in the context of the decorated lattice gas model, 
to a more realistic system, namely a 3D polymer model. Simulation studies of the bond fluc- 
tuation model within the grand canonical ensemble are used to obtain the joint distribution of 
density and energy on the liquid vapour coexistence curve. The critical temperature, chemical 
potential and field mixing parameters of the model are accurately determined by requiring the 
collapse of the measured scaling operator distributions onto their known universal fixed point 
forms. In the sub-critical region close to the critical point, the observed asymmetries of the 
density distribution are again found to be well described by field mixing transformations. 



2 Background 

In this section we provide a brief overview of the principal features of the mixed field FSS 
theory of reference , placing it within the context of the present work. 

The systems we consider are assumed to be contained in a volume L d (with d = 3 in the 
simulations described below) and thermodynamically open so that the particle number can 
fluctuate. The observables on which we shall focus are the particle number density: 

p = L- d N (2.1) 

and the dimensionless energy density: 

u = L- d iw _1 $({r}) (2.2) 

where $({r}) is the configurational energy of the system which we assume takes the general 
two-body form: 

$({>}) = E<Kl r i-rjl)> (2.3) 

with (j)(r) the two-body interaction potential (e.g. square-well or Lennard- Jones) whose asso- 
ciated coupling strength (well-depth) we denote w. 
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Within the formalism of the grand canonical ensemble (GCE), the joint distribution of 
density and energy fluctuations, pi(p,u), is controlled by the reduced chemical potential p 
and the coupling strength w (both in units of ksT). The critical point of the system is located 
by critical values of the chemical potential p c and coupling w c . Deviations of w and p from 
their critical values control the sizes of the two relevant scaling field that characterise the 



critical behaviour [14]. In the absence of the special symmetry prevailing in the Ising model, 



the relevant scaling fields comprise (asymptotically) linear combinations of the coupling and 



chemical potential difference [15 



t = w c — w + s(p — p c ) h = p — p c + r{w c — w) (2.4) 

where r is the thermal scaling field and h is the ordering scaling field. The parameters s 
and r are system-specific quantities controlling the degree of field mixing. In particular r is 
identifiable as the limiting critical gradient of the coexistence curve in the space of p and w. 
The role of s is somewhat less tangible; it controls the degree to which the chemical potential 
features in the thermal scaling field, manifest in the widely observed critical singularity of the 



coexistence curve diameter of fluids E|. 



Conjugate to the two relevant scaling fields are scaling operators Ai and £ , which comprise 
linear combinations of the particle density and energy density ]10|, 11] : 



M = j^gf [p - su] £ = y^st ( n ~~ r P\ ( 2 -5) 

The operator Ai (which is conjugate to the ordering field h) is termed the ordering operator, 
while £ (conjugate to the thermal field) is termed the energy- like operator. In the special case 
of models of the Ising symmetry, (for which s = r = 0), M. is simply the magnetisation while 
£ is the energy density. 

The joint distribution of density and energy is simply related to the joint distribution of 
mixed operators: 

Pl(p, u) = jzipfPLi-M-, £) (2.6) 
Near criticality, and in the limit of large system size, pl(A4,£) is expected to be describable 



by a finite-size-scaling relation of the form [11]: 



where 
and 



p L (M,£) ~ A^A^(A+ 5M , A%6£, A M h, A £ r) (2.7a) 
A e = aeL 1 ^ A M = a M L d -^ A M A + M = A £ A+ = L d (2.7b) 



5M = M- <M> C 5£ = £- <£ > c (2.7c) 
The subscripts c in equations |2.7c| signify that the averages are to be taken at criticality. 



Given appropriate choices for the non-universal scale factors clm an d as (equation |2.7b| ), the 
function p* is expected to be universal. 

Precisely at criticality, equation |2.7a| implies simply 



PL (M,£) ~ A+ A+j5*(A+ 6M,A+8£) (2.8) 

where p*(x, y) = p(x, y, 0, 0) is a function describing the universal and statistically scale invari- 
ant operator fluctuations characteristic of the critical point. 
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In what follows, we shall employ Monte Carlo simulations to explicitly test the prediction 
of equation |2i] for a decorated lattice gas model and a polymer system, both members of the 
Ising universality class. To do so, however, we first require an independent estimate of the 
fixed point function p*(x,y) appropriate to the 3D Ising universality class. In practice, this 
is most readily obtained by considering the prototype member of the Ising class, namely the 
3D Ising model itself. Owing to its lack of field mixing, the scaling operators of the Ising 
model are simply M. — > m (the magnetisation) and £ — > u (the energy density). Moreover, the 
availability of highly accurate estimates for the Ising model critical temperature, circumvents 
the need to perform a time consuming search for the critical point. 



3 Monte Carlo simulations 

3.1 The 3D Ising model and the form of p*(x,y) 

Using a vectorised algorithm on a Cray YMP supercomputer, we have performed high pre- 
cision Monte Carlo simulation measurements of the joint magnetisation and energy density 
distribution for the 3D Ising model on a periodic lattice of side L = 20. The measurements 
were performed at the estimated (reduced) critical coupling K\ = 0.2216595(26), as obtained 
in a previous high precision Monte Carlo study [[|. Following an initial equilibration period 
of 2 x 10 6 Monte Carlo steps per spin (MCS), the magnetisation and energy density were 
sampled at intervals of 50 MCS (in order to reduce correlations), and the results accumulated 
in a histogram. The final histogram (comprising 2 x 10 7 entries) was formed from 12 inde- 
pendent runs, thereby allowing the statistical independence of the data to be assessed and 
statistical errors to be assigned to the results. The resulting form of p*(x,y), normalised to 
unit integrated weight and scaled to unit variance along both axes, is shown in figure [TJ The 
associated ordering operator distribution p*(x) = J p*(x,y)dy , and energy operator distribu- 
tion p*(y) = J p*(x, y)dx are shown in figure [|. One observes that while the form of p*(x) is 
doubly peaked and symmetric, that of p*(y) is singly peaked and asymmetric. 



3.2 The 3D decorated lattice gas model 

The decorated lattice gas model was first proposed in its two dimensional form by Mermin 
as an example of a system exhibiting a singular coexistence diameter [ 17|| . The model was 



subsequently generalised to simple, body, and face centred cubic lattices by other workers 
[ |18| , [19] , |20f1 and studied for its interesting coexistence properties. The simple cubic form of the 
model, on which we shall focus in the present work, consists of an ordinary simple cubic lattice 
gas (whose sites we term the primary sites), augmented (decorated) by additional secondary 
sites on the intersitial bonds. Particles on primary sites are assumed to interact with one 
another via a dimensionless coupling of strength A, while particles on the secondary sites 
interact with those on primary sites via a dimensionless coupling rj, but do not interact with 
each other. A schematic representation of a unit cell of the model is shown in figure |3|. 
The configurational energy <3>({cx}) of the decorated lattice gas model is given by 

$ (W) = VViVj + Aa m cr n (3.1) 

<i,j> [m,n] 

with (Tj = 0, 1. The site indices i and j are taken to run over nearest neighbour primary and 
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secondary sites of the model, while m and n run only over nearest neighbour primary sites. It is 
straightforward to show that the particle-hole symmetry that obtains in the ordinary lattice gas 
is equivalent to the requirement that all sites have the same average energy environment. The 
presence of two inequivalent sublattices in the decorated model clearly violates this condition 
and leads to field mixing. 

Aside from its field mixing properties, the chief asset of the decorated lattice gas model, 
is its analytic tractability. The grand partition function of the asymmetric model can be 
related by means of analytic transforms to that of the ordinary lattice gas model which is itself 
isomorphic to the Ising model. Specifically, one finds: 

T) = (1 + e^ /feT ) 3l7 n(7Z, T) (3.2) 

where Q is the partition function of the decorated model and fi and T are the chemical potential 
and temperature respectively. Bars denote quantities in the ordinary lattice gas and N is the 
number of primary sites in the model. 

Introducing the dimensionless chemical potential £ = [i/k b T equation [J]2] leads to the 
following relationships [19] : 



A 



£ + 61n 
A + ln 



1 + e 



l + e« 
l + e*)(l + e*+ 2 ") 
(1 + e^) 2 



(3.3a) 
(3.3b) 



where £ = ~p/kbT and A = X/kbT are respectively the dimensionless chemical potential and 
dimensionless nearest neighbour coupling constant of the ordinary lattice gas. 

In the 3D ordinary lattice gas the liquid-vapour coexistence line is specified by the condition 
£ = — 3A The location of critical point that terminates this line is not known exactly, but 
is trivially related to that of the Ising model: 



X C = 4K J C , £ C = -3A C (3.4) 

where 

K[ = 0.2216595(26) (3.5) 

is the estimated dimensionless critical coupling of the simple cubic Ising model ||. Estimates 
for the critical point parameters are thus obtainable by feeding this value of K c into equa- 
tions |3.3a| and |3.3b| . Similarly, the coexistence curve of the decorated model is obtainable by 
setting £ = — 3A in equations |3.3a| and |3.3b| to find: 



£ + 3A + 31n 







(3.6) 



l + e« 

Note however, that since A and r\ both enter only as multiplicative factors in the configurational 
energy (equation |3.1| ), the coexistence curve is uniquely parameterised by the value of the 
coupling ratio A/77. Varying this ratio allows one to tune the degree of field mixing [|16| . 
Indeed for the special choice A/77 = —1/3, the average energy environment is identical for 
atoms on both sublattices of the model and particle-hole symmetry is restored. 

Knowledge of the mapping between the decorated lattice gas and the ordinary lattice gas 
also permits an analytic calculation of the field mixing parameters r and s for the model. The 
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value of r is obtainable from equation |3]6| simply by calculating the limiting critical gradient 
of the coexistence curve. A calculation of s proceeds from the observation that in the ordinary 
lattice gas, the field-like scaling field h coincides with the line A = A c in the space of £ and 
A. It follows that in the decorated lattice gas model, the direction of h can be obtained from 
equation |3.3b| by setting A = A c , 77 = 1](X) and solving for A. The value of s is then given by 



(3.7) 



-(I), 

where the derivative is to be evaluated at criticality. 

Compared to more realistic fluid models such as the Lennard- Jones fluid, the great simplic- 
ity of the decorated lattice gas model renders it highly computationally tractable. Moreover 
the prior availability of accurate values for the critical parameters obviates the need for a time 
consuming search of parameter space for the critical point, and simplifies the task of data 
analysis. The model therefore provides an ideal test-bed for simulation studies of critical point 
field mixing. 

3.2.1 The critical limit 

Using a Metropolis algorithm within the grand canonical ensemble (GCE) ||22|| , we have per- 
formed detailed simulation measurements of the joint density and energy distribution p L (p,u) 
of the decorated lattice gas model, at the estimated critical parameters obtained as described 
above. All simulations were performed for a choice of the coupling ratio A/77 = 0.1, for which 
it is known that the coexistence curve of the model closely resembles those of many real fluids 



19fl . In the course of the simulations, three system sizes were studied having linear extent 
L = 12, L = 20 and L = 32. Periodic boundary conditions were employed throughout. Prior 
to data collection, equilibration periods of 5 x 10 6 MCS were utilised. Samples of the density 
and energy density were then performed at intervals of 50 MCS (to reduce correlations) and 
the data stored in histograms. For each system size, 12 independent runs were performed in 
order to test the statistical independence of the data and to assign statistical errors to the 
results. The final histograms of Pl{p, u) comprised 2 x 10 7 entries for the L — 12 and L = 20 
system sizes, and 1 x 10 7 entries for the L = 32 system size. 

The measured form of Pl(p, u) for the decorated lattice gas model is presented in figure £| for 
the L = 20 system size. Clearly apart from a general overall double peaked structure, the form 
of this distribution bears little resemblance to that of p*(x,y) (c.f. figure |I]) which represents 
the joint critical order parameter and energy density distribution in the absence of field mixing. 
To illustrate the differences it is instructive to compare the fluid density distributions Pl(p) = 
J pi{p,u)du with p*(x), and the fluid energy density distribution Pl(u) = J pL(p,u)dp with 
p*(y) . 

The forms of Pl{p) for all three system sizes are shown in figure |5|a. In contrast to p*(x) 
(figure ^|a), all the density distributions exhibit a pronounced asymmetry, qualitatively similar 
in form to that observed in the 2D version of the same model |I6| . Even more conspicuous, 



however, are the differences between the finite-size forms of Pl{u) (figure |5|b), and that oip*(y) 
(figure Qb). Clearly while the latter is singly peaked, the former are doubly peaked, with a 
form (were one to plot Pl{~u))., reminiscent of the density distribution. As we shall show, the 
explanation of these differences is to be found in the field mixing that manifests the lack of 
particle-hole symmetry in fluids. 
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In order to expose the universality linking the critical point of the decorated lattice gas 
model to that of the Ising model, it is necessary to recast pl(p,u) in terms of the joint dis- 
tribution of scaling operators pl(M.,£) , cf. equation |2~8] . To do so however, requires specifi- 
cation of the field mixing parameters s and r featuring in the definitions of M. and £ (equa- 



tion [2.5|) . In practice, these values may be readily found by requiring that the single scaling 
operators distributions Pl(-M) and Pl{£) match their respective fixed point forms p*(x) and 
p*(y) ■ Carrying out this procedure yields the matchings shown in figure |||a and |||b. The 
associated estimates for the field mixing parameters are s = —0.143(8) and r = —3.11(1). 
These values compare very favourably with those calculable analytically, for which one finds 
s = —0.1428 • • • , r = —3.1163 • • • Such a high level of accord indicates that the matching of 
the operator distributions to the universal Ising forms is a potentially very accurate method 
for determining the field mixing parameters. 

Having obtained estimates for s and r, one may then construct the joint distribution of 
scaling operators pl(A4,£). The resulting form is shown in figure 0, and should be compared 
with that of p*(x,y) shown in figure |1|. Clearly the agreement between the operator distribu- 
tions and the universal fixed point form is gratifying, providing substantial corroboration of 
the mixed field FSS theory. 

3.2.2 The subcritical region 

As the critical point is approached along the line of phase coexistence, the known symmetries 
of the Ising problem imply that 

(M)* - (M) c = ±a\rf, (3.8a) 
(£}^ — (£) c = &|r| 1_a + terms analytic at criticality (3.8b) 

where a and b are critical amplitudes and ± denote limits as the coexistence curve is ap- 
proached from above (h — > + ) or below (h — > 0~). Recalling that M. = p + s£ then yields the 
two branches of the coexistence curve densities near criticality: 

p± — p c = ±a\rf + s6|r| 1_a + terms analytic at criticality, (3.9) 
which displays a singular diameter: 

Pd~Pc= \{P+ + P-) ~Pc~ Shir] 1 -* (3.10) 



as is indeed observed experimentally ||13[| . 

At temperatures outside the critical region, the relations |3.8a| and |3.8b| are not generally 



expected to hold and a crossover description to regular classical behaviour is more appropriate 



^|, |27|, ^8) . Nevertheless as observed in many computer simulations of various simple 
fluids, the temperature dependence of the order parameter is quite accurately described by Ising 
critical exponents over a remarkably wide range of subcritical temperatures |29|| , without the 
need to introduce higher order (non-linear) terms in the Wegner expansion of the scaling fields 
TM. In view of this it seems of interest to examine the range of applicability of the scaling 



form |2.7a| in the subcritical two-phase region. 

To this end we have obtained the joint density and energy density distribution of the deco- 
rated lattice gas model at temperatures T = 0.9T C and T = 0.8T C , for a system of side L = 20. 
In order to circumvent the prohibitively large 'tunneling' times between the coexisting phases 
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that normally plague GCE simulations in the two-phase region, we have employed the multi- 



canonical preweighting scheme j30|. This scheme uses weighted transitions to encourage the 
system to sample those interfacial configurations that would otherwise occur only very rarely. 
The weight factors are chosen such that the sampled density distribution is approximately flat 
in the density region between the peaks of the two coexisting phases. After the simulation, the 
correct ('unweighted') coexistence distribution is regained by dividing out the weight factors 
from the sampled distribution. In this manner the effective tunneling frequency between the 
coexisting phases may be increased by many orders of magnitude, thereby facilitating very 
accurate estimates of the coexistence form of pi(p,u). 

Using the multicanonical preweighting scheme, samples of the density and energy were ac- 
cumulated every 20 MCS and stored in histograms. The final histograms for Pl{p, u) comprised 
approximately 5 x 10 7 entries. In figure |S]a we present the measured coexistence forms of Pl{p) 
for T = 0.9T C and T = 0.8T C . Also included in the figure is the critical point form of p^(p) for 
the L = 20 system size, as obtained previously (c.f. subsection |j.2.1p . Clearly the distributions 
are all asymmetric, those for the two subcritical temperatures having a form similar to those 
observed in an asymmetric version of the 2D Blume-Emery-Griffiths spin model fl3l| . Although 
the two peaks of the sub-critical distributions have equal weight, (reflecting the thermodynamic 
condition for coexistence ||32||), one observes that the high density peak is much broader and 
shorter than the low density peak. Moreover the magnitude of the change in position of the 
high density peak on lowering the temperature is much greater than that in the low density 
peak. It is this effect that gives rise to the asymmetric form of the temperature-density phase 
diagram of fluids. 

In figure |8|b we present the corresponding forms of the ordering operator distribution 
Ph{Ai). For each temperature the value of the field mixing parameter r was obtained as 
the gradient of the coexistence curve, while the value of s was obtained by applying the pre- 
scription of equation |3]7]. One sees that in this choice of variables, the distributions at T = T c 
and T = 0.9T C are largely symmetric i.e. both peaks have the same height and width. Only the 
distribution for T = 0.8T C exhibits small deviations from symmetry, these being attributable to 
the non-singular (and non-Ising) part of the partition function (equation 3^). Thus it would 
appear that in the present case at least, the validity of the scaling form |2.7a| extends some 
10% or more below the critical temperature. In this temperature range, the asymmetry of the 
density distribution can therefore be accurately ascribed to the linear mixing of the energy 
density into the ordering operator. It remains to be seen however, to what extent this finding 
holds true in more realistic fluid models as well as those possessing very large coexistence curve 
asymmetries, such as polar |34| , ionic [[£J or metallic fluids [3^ 



It is also instructive to compare the simulation results with analytic calculations of the 
coexistence density diameter pd and the ordering operator diameter Aid- The latter can 
be calculated exactly, while the former may be obtained approximately by employing the 
tabulated values of the Pade Approximants for the temperature dependence of the Ising model 
energy density |33|]. Performing these calculations for fractional temperatures t = T/T c in 
the range 0.7 < t < 1.0 yields the results shown in figure |^. The weak critical singularity in 
the coexistence diameter is readily discernible from the figure although, as expected, no such 
singularity obtains for Aid{t), which is analytic at the critical point. Some variation in Aid 
is seen as a function of temperature, however, but this again arises from the non-singular, 



non-Ising part of the partition function |3.2| . Also included in the figure are the simulation 



estimates of (Ai)(t) and (p)(t) obtained from the multicanonical simulations at T = 0.9T C and 
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T = 0.8Tc, as well as the conventional simulations at the critical point. Clearly a very good 
overall agreement between the simulations and the analytical predictions is apparent. 



3.3 The universal critical finite-size spectrum of 

p L (p) and p L (u) 

In this subsection we return to a consideration of the critical point forms of Pl{p) and Pl{v) 
in fluids, with the aim of gaining an understanding of their shapes and finite-size-scaling 
behaviour. To this end it is expedient to reexpress p and u in terms of the scaling operators. 



Appealing to equation 2.5, one finds 



u = 8-rM p = M-sS, (3.11) 
so that the critical density and energy density distributions are 

p L (u) = p L (£ - rM) p L (p) = Pl(M - sS) (3.12) 

Now the structure of the scaling form |2.7aj shows that the typical size of the fluctuations in 
the energy-like operator will vary with system size like 58 ~ [ J -( 1 - a )/ u ) while the typical size 
of the fluctuations in the ordering operator vary like 5M ~ L~^ v . It follows that for a given 
L, the shape of the energy and density distributions can be identified with the distribution of 
the variable 

X e = aj}5M cosG + a^S£ sin 6, (3.13) 

with 

tan0u = _^_2T{i-«-fl/" andtan0 p = ^ L -(^-^ ( 3 .14) 
ra M a M 

where the subscripts u and p signify that the value of corresponds to the energy density and 
density distributions respectively. 

The distributions p(X&) constitute a one-parameter class of universal functions describing 
the density and energy distributions of fluids at finite L. Geometrically, can be interpreted as 
defining a direction OX@ in the basal plane formed by the Ox and Oy axes of figure [l], making 
an angle with the Ox axis. The form of p(X®) is then obtainable by projecting p*(x, y), onto 
the vertical plane which includes the line OX®. A representative selection of such projections 
is shown in figure ITUI. For = 0° one obtains simply the ordering operator distribution p*(x) , 



while for = 90° the from is that of the energy-like operator p*(y) distribution. Intermediate 
between these values a range of behaviour is obtained, representing the finite L forms of Pl{p) 
and Pl(u)- 

Asymptotically (i.e. as L — > oo), equation [3.14| implies that both U and Q p approach zero 



4(| so that in this limit 

Pl(u) = p L (-rM) ~ aM-W^p^-aM-'rL^SM) (3.15) 

Pl(p) = p L (M) ~ a M - l LV»p* M (a M - l LP/»5M) (3.16) 

It follows that for any finite s and r, the limiting critical point forms of Pl(p) and Pl{u) 
both match the critical ordering operator distribution p*(x). The approach to this limiting 
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behaviour is indeed clearly evident in the distributions of figure |5]. We note however that the 
limiting form of Pl{u) differs radically from that of the Ising model where, owing to the absence 
of field mixing, lini^oo Pl{u) = P*(y)- The profound influence of field mixing on the critical 
behaviour of fluids should therefore be apparent [[41] . 



Finally in this subsection, we point out that precise knowledge of the location of the critical 
point (obtained eg. from the data collapse of the scaling operators onto their fixed point 
forms), does not imply the possibility of directly measuring the infinite- volume density and 
energy density. To appreciate this, recall that 

(u) c = (£) c - r(M) c (p) e = (M) e - s(£) c (3.17) 

Now, while symmetry considerations dictate that the value of (A4) c = J Pl{M)<IM. is inde- 
pendent of system size, no such symmetry condition pertains to Pl{£), whose average value 
(£) c = J d£p L (£) at criticality is expected to vary with system size like 

(£) C (L) - (£> c (oo) ~ (3.18) 

It follows that in order to extract infinite volume estimates of p c and u c from simulations at 
the critical point, it is necessary to extrapolate data from a number of different system sizes 



to the thermodynamic limit. This procedure is illustrated in figure [TT]for the decorated lattice 
gas model, using critical point data from the three system sizes L = 12, 20, 32. The measured 
values of (p) c (L) and (u) c (L) are plotted against L~( 1 ~ a ^ u . A least squares fit to the data 
yields the infinite volume estimates p c = 0.3371(1) and u c = —0.8385(6). We remark that 



the existence of these finite-size shifts imply that the equal weight criterion |32|, |31] for the 
order parameter distribution, while correctly identifying the coexistence curve in the subcritical 
regime, must fail close to the critical point pT 



3.4 Liquid-vapour equilibria of a polymer model 

We now apply the techniques developed in the foregoing sections to the study of critical 
phenomena and phase coexistence in a more realistic model fluid, namely a polymer system. 
The model we consider is the bond fluctuation model, a coarse grained lattice-based polymer 
model which combines the essential qualitative features of real polymer systems — monomer 
excluded volume and connectivity — with computational tractability. Within the framework of 
the model, each monomer occupies a whole unit cell of a periodic simple cubic lattice. Excluded 
volume interactions are catered for by requiring that no lattice site can simultaneously be 
occupied by two monomer corners. Monomers along the polymer chains are connected by 
bond vectors which can assume one of 108 possible values, providing for 87 distinct bond 
angles and 5 distinct bond lengths. For a more detailed description of the model, the reader 
is referred to the literature |Tj] . 

Using a grand canonical simulation algorithm, we have simulated chains of length N = 20 
monomers, interacting via a short range square well potential, the range of which was set at a/6 
(in units of the lattice spacing). Chain insertions and deletions were facilitated by use of the 
configurational bias Monte Carlo (CBMC) method of Siepmann |43 |. The essential idea behind 



the CBMC method is to improve the low acceptance rate associated with random trial chain 
insertion, by 'growing' chains of favourable energy into the system. A bookkeeping scheme 
maintains a record of the statistical bias associated with choosing favourable chain conforma- 
tions, and this bias is subsequently removed when the acceptance probability is calculated. 
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The CBMC technique has also recently been used in conjunction with Gibbs ensemble Monte 



Carlo simulations of liquid- vapour phase coexistence of off-lattice alkanes models [44 
The quantities measured in the simulations were the monomer density: 

p = 8nN/V (3.19) 

and the dimensionless energy density: 

u = 8^" 1 $({r})/K (3.20) 

where n is the number of chains, $({r}) is the configurational energy, w is the well depth and 
V is the system volume which was set at V — 40 3 . Here the factor of 8 derives from the number 
of lattice sites occupied by one monomer. In the course of the simulations, measurements of 
p and u were performed at intervals of 5000 chain insertion attempts and accumulated in the 
joint histogram pl(p,u). The final histogram comprised some 3 x 10 5 entries. 

In contrast to the decorated lattice gas model considered in previous sections, the line 
of liquid-vapour phase coexistence is not known a-priori for the polymer system and must 
therefore be identified empirically. The precise location of the coexistence curve is prescribed 
by the equal weight criterion for the two peaks of the density distribution |32| . Unfortunately, 
the task of identifying the coexistence curve using this criterion is an extremely time consuming 
and computationally demanding one, since the density distribution is generally very sensitive 
to small deviations from coexistence. In practice, however, it suffices to obtain data for only a 
few points close to the coexistence curve. The full coexistence curve between these points can 



subsequently be constructed using histogram reweighting techniques |45], [4(|. Provided that 
the measured density distributions are doubly peaked and the temperatures studied are not 
too widely separated, this technique permits a very accurate determination of the coexistence 
curve locus. 

Starting with an initial well-depth w = 0.569, the approximate value of the coexistence 
chemical potential was determined by tuning p until the density distribution exhibited two 
peaks. Again the multicanonical preweighting scheme |3(J was employed in order to overcome 
the otherwise very large tunnelling times between the coexisting phases. A histogram extrap- 
olation based on this data was then used to estimate the value of the coexistence chemical 
potential for a well-depth w = 0.56, which lies close to the critical well-depth. A further 
long runs was carried out at this near-coexistence point. By extrapolating the measured near- 
coexistence histograms of pl(p,u) in conjunction with the equal weight criterion, we were 
then able to construct a sizeable portion of the coexistence curve (and its analytic extension). 
Representative forms of the density distributions along the line of coexistence and its analytic 
continuation are shown in figure |l2|a. The coexistence curve, expressed as a function of the 
well-depth w and chemical potential p is shown in figure [L3|. 

To locate the critical point along the line of phase coexistence, we utilised the universal 
matching condition for the operator distributions Pl{M) and Pl{£) ■ Again applying the his- 
togram reweighting technique, the well-depth, chemical potential and field mixing parameters 
were tuned until the forms of pl(M) and Pl{£) most accurately matched the universal critical 
Ising forms of figure |2|. The results of performing this procedure are shown in figure |14[ Given 
that the system contains an average of only about 100 polymer chains, the quality of the data 
collapse is remarkable. The mappings shown were effected for a choice of the parameters 

w c = 0.5584(1), p c = -5.16425(2), s = -0.135(4), r = -2.55(2) (3.21) 
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where we have defined p, to be the chemical potential per monomer. The corresponding critical 
density and energy density distributions are shown in figure |15[ They yield the (finite-size) 
estimates p c = 0.199(3) and u c = -0.304(4). 

Turning lastly to the subcritical two-phase regime, we have again considered the effect of 
forming linear combinations of the density and energy density. Figure |T4]b shows the form of the 
ordering operator distribution Pl{M) obtained at coexistence using the histogram reweighting 
technique for the same values of w shown in figure 14 a. In each instance the values of s was 
chosen so that the two peaks of pl(M) had both equal heights and equal weights, while the 
value of r was chosen so that Pl(£) was singly peaked. As was the case for the decorated 
lattice gas model, simple field mixing transformations also appear to account for the sub- 
critical coexistence curve asymmetries of the polymer density distribution, at least over the 
limited range of w studied here. 



4 Conclusions 

In summary we have provided explicit demonstration of the field mixing transformations that 
link the fluctuation spectra of the order parameter and energy in the critical fluid to those of 
the critical Ising magnet. The results serve to underline the profound influence of field mixing 
on the non universal critical behaviour of fluids. This influence is manifest most notably as a 
finite-size shift to the measured critical density, and as an alteration to the limiting (large L) 
form of the critical energy distribution. Field mixing is also found to account for the observed 
asymmetries of the coexistence density distribution over a sizeable portion of the sub-critical 
region. 

With regard to the general computational issues raised in this study, it has been seen that 
effecting the data collapse of the fluid scaling operator distributions onto their (independently 
known) universal fixed point forms, provides a very powerful method for accurately locating 
the critical point and determining the field mixing parameters of model fluids. This use of the 
scaling operator distributions represents the natural extension to fluids of the order parameter 
distribution method deployed so successfully in the study of symmetric spin models. Thus, in 
principle at least, there would appear to be no barriers to attaining similar degrees of accuracy 
in the study of critical fluids as has previously been achieved for lattice spin systems. 

The successes of the present work (and of an earlier FSS study of the 2D Lennard- Jones 
fluid |l IJ ) also attest to the utility of the grand canonical ensemble for simulation studies of 
near-critical fluids. The benefits of this ensemble stem principally from the fact that density 
fluctuations are observable on the scale of the system size itself, thus freeing the method of 
the interfacial effects and additional length scales that complicate use of the 'sub-block' finite- 
size scaling technique within the canonical (NVT) ensemble || |47| . High quality results can 
therefore be obtained using comparatively much smaller system sizes, with concomitant savings 
in computational effort. 

The ability to perform a full finite size scaling analyses in the near-critical region also rep- 
resents an important advantage of the GCE approach over the Gibbs ensemble Monte Carlo 
(GEMC) simulation technique f29j. In the GEMC method, the fluctuating box size seems to 
preclude a rigorous FSS analysis ||S, 4"9|] , thus seriously hindering the accurate location of the 



critical point. The GEMC is, nevertheless, very efficient for locating the temperature-density 
phase diagram in the sub-critical regime. For this task, use of the bare GCE method is only 
feasible for temperatures within a few percent of the critical temperature because the other- 
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wise high interfacial free energy results in prohibitively large 'tunneling' times between the 
coexisting phases. Nonetheless, as we have shown, this problem is surmountable by combining 
the GCE with recently developed multicanonical preweighting and histogram reweighting tech- 
niques, thereby enabling accurate studies of the coexistence density and energy fluctuations 
even well below the critical temperature. 

Acknowledgements 

The authors have benefitted from helpful discussions with B.A. Berg, K. Binder and A.M. 
Ferrenberg. NBW acknowledges the financial support of a Max Planck fellowship from the 
Max Planck Institut fiir Polymerforschung, Mainz. Part of the simulations described here 
were performed on the CRAY-YMP computers at the HLRZ Jiilich and the RHRK Universitat 
Kaiserslautern. Partial support from the Deutsche Forschungsgemeinschaft (DFG) under grant 
number Bi3 14/3-2 is also gratefully acknowledged. 

References 

[1] For a review, see V. Privman (ed.) Finite size scaling and numerical simulation of statis- 
tical systems (World Scientific, Singapore) (1990). 

[2] K. Binder, Z. Phys. B 43, 119 (1981). 

[3] A.M. Ferrenberg, D.P Landau, Phys. Rev. B 44, 5081 (1991) 

[4] D. Nicolaides and A.D. Bruce, J. Phys. A. 21, 233 (1988). 

[5] A.P Gottlob and M. Hasenbusch, Physica A201, 593 (1993). 

[6] P. Peczak, A.M. Ferrenberg and D.P. Landau, Phys. Rev. B43, 6087 (1991). 

[7] K. Chen, A.M. Ferrenberg and D.P. Landau, Phys. Rev. B48, 3249 (1993). 

[8] K. Binder in Computational Methods in Field Theory H. Gausterer, C.B. Lang (eds.) 
Springer- Verlag Berlin-Heidelberg 59-125 (1992). 

[9] M. Rovere, D.W. Heermann and K. Binder, J. Phys. Condens. Matter 2, 7009 (1990). 

[10] A.D. Bruce and N.B. Wilding, Phys. Rev. Lett. 68, 193 (1992). 

[11] N.B. Wilding and A.D. Bruce, J. Phys. Condens. Matter 4, 3087 (1992). 

[12] Following convention we employ the language of the ordinary lattice gas model, which is 
formally equivalent to the Ising model. 

[13] J.V. Sengers and J.M.H. Levelt Sengers, Ann. Rev. Phys. Chem. 37, 189 (1986). 

[14] F.J. Wegner, Phys. Rev. B. 5, 4529 (1972). 

[15] J.J. Rehr and N.D. Mermin, Phys. Rev. A 8, 472 (1973). 

[16] N.B. Wilding, Z. Phys. B93, 119 (1993). 



13 



[17] N.D. Mermin, Phys. Rev. Lett. 26, 957 (1971) . 

[18] J.A. Zollweg and G.W. Mulholland, J. Chem. Phys. 57, 1021 (1972). 

[19] G.W. Mulholland and J.J. Rehr, J. Chem. Phys. 60, 1297 (1974) 

[20] G.W. Mulholland, J.A. Zollweg and J.M.H. Levelt-Sengers, J. Chem. Phys. 62, 2535 
(1975). 

[21] T.D. Lee and C.N. Yang, Phys. Rev. 87, 410 (1952). 

[22] K. Binder and D.W. Heermann, Monte- Carlo Simulation in Statistical Physics, 2nd ed. 
(Springer, Heidelberg, 1992). 

[23] R.R. Singh and K.S. Pitzer , J. Chem. Phys. 90, 5742 (1989). 

[24] J. Luettmer-Strathmann, S. Tang and J.V. Sengers, J. Chem. Phys. 97, 2705 (1992). 

[25] M.A. Anisimov, S.B. Kiselev, J.V. Sengers and S. Tang, Physica A188, 487 (1992). 

[26] Z.Y. Chen. PC. Albright and J.V. Sengers, Phys. Rev. A41, 3161 (1990). 

[27] S. Tang, J.V. Sengers and Z.Y. Chen, Physica A179, 344 (1991). 

[28] Z.Y. Chen, A. Abbaci, S. Tang and J.V. Sengers, Phys. Rev. A42, 4470 (1990). 

[29] A.Z. Panagiotopoulos, Mol. Phys. 61 813 (1987); A.Z. Panagiotopoulos, Mol. Simul. 9 1 
(1992). 

[30] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992). 

[31] C. Borgs and S. Kappler, Phys. Lett. A171, 37 (1992); C. Borgs, P.E.L. Rakow and S. 
Kappler, J. Phys I (France) 4, 1027 (1994). 

[32] C. Borgs, R. Kotecky, J. Stat. Phys 60, 79 (1990), Phys. Rev. Lett 68, 1734 (1992), 

[33] PE. Scesney, Phys. Rev. Bl, 2274 (1970). 

[34] M.E. Van Leeuwen, Mol. Phys. 82, 383 (1993). 

[35] J.M.H. Levelt Sengers and J.A. Given, Molec. Phys. 80 899 (1993). 

[36] S. Jiingst, B. Knuth and F. Hensel, Phys. Rev. Lett. 55, 2160 (1985). 

[37] R.E. Goldstein, A. Parola, NW. Ashcroft, M.W. Pestak, M.W.H. Chan J.R. de Bruyn 
and D.A. Balzarini, Phys. Rev. Lett. 58, 41 (1987). 

[38] R.E. Goldstein, and A. Parola, J. Chem. Phys. 88, 7059 (1988). 

[39] M. Miiller and N.B. Wilding, Mainz preprint. 

[40] The universal function representing the difference between the finite-size form of Pl{p) 



and its limiting form p*(x) has been previously studied in references |Tl], [T6[] , where it was 
termed the 'field mixing correction'. 

14 



[41] The differences between the energy distributions of fluids and symmetric systems have 
also recently been discussed by us elsewhere in the context of the critical finite-size scaling 
behaviour of the heat capacity in a binary polymer mixture |3!| . 

[42] I. Carmesin and K. Kremer, Macromolecules 21, 2819 (1988); H-.P. Deutsch and K. 
Binder, J. Chem. Phys. 94, 2294 (1991) 

[43] J.I. Siepmann, Mol. Phys. 70, 1145 (1990) 

[44] J.I. Siepmann, S. Karaborni and B. Smit, J. Am. Chem. Soc. 115, 6454 (1993); Nature 
365, 330 (1993). 

[45] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61 2635 (1988); A.M. Ferrenberg 
and R.H. Swendsen Phys. Rev. Lett. 63, 1195 (1989). 

[46] H.-P. Deutsch, J. Stat. Phys. 67, 1039 (1992). 

[47] M. Rovere, P. Nielaba and K. Binder, Z. Phys. B90, 215 (1993). 

[48] J.R. Recht and A.Z. Panagiotopoulos, Mol. Phys. 80, 843 (1993); A.Z. Panagiotopoulos, 
preprint. 

[49] K.K. Mon and K. Binder, J. Chem. Phys. 96, 6989 (1992). 



15 



Figure 1: Estimates of the fixed point form of the joint scaling operator distribution p*(x,y) 
appropriate to the 3D Ising universality class, obtained as the joint magnetisation and energy 
density distribution of the L = 20 periodic 3D Ising model at the estimated critical coupling 
K = 0.2216595. The data is expressed in terms of the scaling variables x = a~j^L^l v {M. — M. c ) 
and y = a^ 1 L^~ a ^ u (£ — £ c ), with the scale factors a^ 1 and aj^ chosen so that the distribution 
has unit variance along both axes. 

Figure 2: (a) The ordering operator distribution p*(x) appropriate to the 3D Ising universality 
class, obtained as the magnetisation distribution of the L = 20 3D periodic Ising model at the 
assigned value of the critical coupling K* = 0.2216595, and expressed in terms of the scaling 
variable x = ajJ t L l3 ^ l/ (M — M c ). (b) The corresponding energy operator distribution p*(y), 
expressed in terms of the scaling variable y = a E x lS v ~ OL ^ 1 ' [E — S c ). In both cases, the value 
of the non-universal scale factors a^ 1 and aj^ were chosen so that the distributions have unit 
variance. Statistical errors do not exceed the symbol sizes. 



Figure 3: Schematic representation of a unit cell of the decorated lattice gas model. The 
primary sites have been shaded, while the secondary (decoration) site have not. Primary 
sites interact with other nearest neighbour primary sites with a coupling A, and with nearest 
neighbour secondary sites with a coupling r\. 
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Figure 4: Estimates of the normalised joint distribution of density and energy density Pl(p, u) of 
the 3D decorated lattice gas model for L = 20, at the assigned value of the critical parameters. 

Figure 5: (a) The density distributions of the 3D decorated lattice gas model at the assigned 
values of the critical parameters, for system sizes L = 12, L = 20 and L = 32. (b) The 
corresponding energy density distributions. Statistical errors do not exceed the symbol sizes. 

Figure 6: (a) The critical ordering operator distribution Pl{-M) of the decorated lattice gas 
model for the system sizes L = 20 and L = 32. The L = 12 data has been omitted for clarity. 
Also shown for comparison is the universal ordering operator distribution p*(x) (c.f. figure |2|a). 
(b) The critical energy operator distribution Pl{£) compared with the universal fixed point 
form p*(y) (figure 0b). In both cases, the non-universal scale factors a^ 1 and aj^ have been 
chosen to give unit variance for the L = 32 data set. The exponent ratios were taken to be 
(3/v — 0.5176 and ajv = 0.177. Statistical errors do not exceed the system sizes. 

Figure 7: The joint critical distribution of scaling operators pl(A4,£) for the 3D decorated 
lattice gas model, obtained by applying appropriate field mixing transformation to the distribu- 
tion of figure B. The data is expressed in terms of the scaling variables x = aJ^L /3 ^ u (A4 — A4 C ) 
and y = a £ l lJ 1 ^ a ^ p '{£ — £ c ). The values of the non-universal scale factors aj 1 and aj^ have 
been chosen to give unit variance along both axes. 



Figure 8: (a) Estimates of the coexistence density distributions of the L = 20 decorated lattice 
gas model at T = T c , T = 0.9T C and T = 0.8T C . (b) The corresponding forms of the ordering 
operator distribution pl(A4). Statistical errors do not exceed the symbol sizes. 



Figure 9: Comparison of the calculated (lines) and measured (points) temperature dependence 
of the coexistence density and ordering operator diameters. The calculated temperature de- 
pendence was obtained from exact and series expansion methods, as described in the text. The 
data points were obtained from the simulations, except the estimate for p c , which derives from 
the finite-size extrapolation of figure 11. Statistical errors do not exceed the symbol sizes. 



Figure 10: Selections from the universal finite-size spectrum of critical density and energy 
density distributions of fluids. The distributions were obtained according to the procedure 
described in the text. The values of the non-universal scale factors aj 1 and aj^ have been 
chosen to ensure that the distributions have unit variance 



17 



Figure 11: (a) The measured average density (p) c (L) of the decorated lattice gas model at 
the critical point, expressed as a function of L~( 1 ~ a ^ u . The least-squares fit yields an infinite 
volume estimate p c = 0.3371(1). (b) The measured average energy density (u) c (L) of the 
decorated lattice gas model at the critical point, expressed clS db function of L~( 1-a )/" . The 
least-squares fit yields an infinite volume estimate u c = —0.8385(6). 



Figure 12: (a) The monomer density distribution for a selection of well-depths w along the 
line of liquid- vapour coexistence. The distributions were obtained by applying the histogram 
reweighting technique to measured near-coexistence data for w = 0.56 and w = 0.569. (b) 
The corresponding distributions of the ordering operator pl (Ai ) . 



Figure 13: The line of liquid- vapour phase coexistence and its analytic extension in the space of 
/i and w, for values of w in the range 0.555 — 0.575. The results were obtained by implementing 
the equal weight criterion for the density distribution. Also shown are the measured directions 
of the relevant scaling fields. 



Figure 14: (a) The ordering operator distribution Pl{M) of the polymer model at the assigned 
critical parameters. Also shown for comparison is the universal fixed point form p*(x) (cf. 
figure ^|a). (b) The energy-like operator distribution Pl{£) compared with the universal fixed 
point form p*(y) (cf. figure 0b). In both cases the non-universal scale factors aj 1 and have 
been chosen to give unit variance. Statistical errors are indicated by representative error bars. 



Figure 15: (a) Estimates for the polymer monomer density distribution Pl(p) at the assigned 
values of the critical well-depth and chemical potential, (b) Corresponding estimates for the 
polymer energy density distribution Pl{u). Representative error bars are shown. 
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